Self-consistent inhomogeneous steady states in Hamiltonian mean field dynamics 
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Long-lived quasistationary states, associated with stationary stable solutions of the Vlasov equa- 
tion, are found in systems with long-range interactions. Studies of the relaxation time in a model 
of TV globally coupled particles moving on a ring, the Hamiltonian Mean Field model (HMF), have 
shown that it diverges as A''^ for large A'', with 7 ~ 1.7 for some initial conditions with homo- 
geneously distributed particles. We propose a method for identifying exact inhomogeneous steady 
states in the thermodynamic limit, based on analysing models of uncoupled particles moving in an 
external field. For the HMF model, we show numerically that the relaxation time of these states 
diverges with with the exponent 7 ~ 1. The method, applicable to other models with globally 
coupled particles, also allows an exact evaluation of the stability limit of homogeneous steady states. 
In some cases it provides a good approximation for the correspondence between the initial condition 
and the final steady state. 
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I. INTRODUCTION 

Systems with long-range interactions exhibit unusual 
thermodynamic and dynamical properties [1] . These sys- 
tems are characterized by two-body interaction poten- 
tials which decay at large distances r as r~", with a 
smaller than the spatial dimension d. With such po- 
tentials these systems are not additive and, as a result, 
they display characteristic features such as ensemble in- 
equivalence, negative specific heat, temperature jumps 
and broken ergodicity. 

Studies of the dynamics of a class of models have shown 
that some initial macroscopic states display slow relax- 
ation to equilibrium. Starting from these initial states the 
system evolves on a short timescale towards a state that 
is not the one predicted by equilibrium statistical me- 
chanics. The lifetime of this quasistationary state (QSS) 
diverges algebraically with system size. For a finite sys- 
tem one eventually relaxes towards thermal equilibrium 
0. The full nonlinear dynamical relaxation process re- 
mains to be fully characterized. 

A prototypical model for which this dynamical behav- 
ior has been explored in detail is the Hamiltonian Mean 
Field model (HMF)j3i], which describes the motion on a 
ring of N fully coupled particles. Its Hamiltonian is 
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where — tt < 9i < tt is the position of the i-th particle 
and Pi is its conjugated momentum. The model has a 
second order phase transition at the energy Uc = H/N ~ 
3/4 at which the order parameter m = {nix, my) — 



J2i{cosdi,sm6i)/N changes from zero (above Uc) to non 
zero (below Uc). This order parameter is referred to as 
the magnetization vector. 

For this mean-field model, as for other long-range inter- 
acting systems, the time evolution of the single-particle 
distribution function f{9,p,t) is conveniently studied in 
the large N limit by the Vlasov equation Q 
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is the effective potential seen by the particles. 

It has been proposed that linearly stable stationary so- 
lutions of the Vlasov equation are associated with QSS of 
the large but finite N system [2] . This has been discussed 
in detail for homogeneous, 9 independent, steady states 
where m = . Studying inhomogeneous (m 7^ 0) steady 
states in this class of models is less straightforward 
Exploring these stationary states, their stability and the 
relaxation process for quasistationarity is a challenging 
goal. 

The stationary solutions fs{0,p) of the Vlasov equa- 
tion are functions of the integrals of motion and thus, in 
one dimension, they are generically of the form 



fs{9,p) = F{h{9,p)) with h{9,p)^ 



p- 
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Here h{9,p) is the single-particle energy. The function 
F has to satisfy a self-consistency condition obtained by 
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inserting Q into ([3]). However, this is satisfied by a wide 
class of functions F, reflecting the fact that the Vlasov 
equation has many stationary states. In this regard, a 
widely used concept is that of BGK modes [1, @|. An 
application of this method to the HMF model has been 
proposed in f7|. We will here propose a different ap- 
proach. 

In this paper we introduce a simple method for con- 
structing inhomogeneous steady states by analyzing the 
dynamics of models of non-interacting particles moving 
in an external potential. Relaxation in such systems 
has already been considered in Refs. however it 

has not been used to derive analytically inhomogeneous 
steady states of the interacting particles systems. A 
different approach yielding approximate inhomogeneous 
steady states in systems with long-range interactions has 
previously been introduced by Lynden-Bell [ll[ and sub- 
sequently applied in a variety of models [lol. ll2| - [l^ | . Our 
method yields exact steady states and enables one to 
evaluate their stability limit. It also provides a procedure 
to relate in an approximate way an initial distribution to 
the steady state to which it evolves. For the sake of clar- 
ity, the method is applied to the HMF model, although 
it could be readily formulated in a more general context. 

The paper is organized as follows. In Section [H] we in- 
troduce the uncoupled particles model and we derive the 
invariant measure. Section Hill presents the application of 
our method to the HMF model. In Section ITVl we derive 
the exact stability limit of the homogeneous waterbag 
state. Section |V] is devoted to some conclusions. 



II. THE UNCOUPLED PARTICLES MODEL 

In order to find the stationary states of the HMF 
model, we consider the dynamics of an ensemble of un- 
coupled particles moving in a fixed external field H point- 
ing, say, in the x-direction. The energy of a single particle 
is given by 



e{9,p) = y -HcosO 



(5) 



For an arbitrary function F{£{9,p)) to be a steady state 
of the interacting model ([1]), H has to satisfy the self- 
consistency condition 



H = rrir = 
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In the following we suppress the subscript x. Any F that 
satisfies these conditions yields an exact stationary so- 
lution of the HMF model. In order to relate an initial 
distribution to the steady state to which it evolves, we 
consider an initial distribution of particles fo{d,p)- The 
dynamics of the uncoupled particles model ([5]) is such 
that particles in a given energy shell [e, e -I- de] keep mov- 
ing inside that shell, eventually reaching a homogeneous 
distribution within it. As a result, the system attains the 



following steady state distribution 
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For simplicity, we present below the analysis for the often 
studied waterbag initial condition 
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In order to evaluate P{9,p), it is convenient to first con- 
sider the energy distribution Pei^)- For the waterbag 
initial state ([8]) it is given by 



-1 /• r^P 2 



Carrying out the integral over p, one gets 
1 /■ .„ 1 



d9 
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ioT —H < e < Ap? 12 — HcosA9 and zero outside this 
range. The integration over 9 needs to be done in the do- 
main enclosed by the waterbag initial condition, namely 



Thus, 



< e + Hcos9 < 
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where 6*0 and 9i satisfy (see Fig. [T^) 



_ arccos(-e/iJ) , for -H<e<-H cos A9 
~ A9 , for e > -iJcosAe* , 

(13) 

, for - iJ < e < Ap^/2 - H 

^ ^ J arccos C^^y-' ) , for Ap^/2 -H<e< 
^ ^ Ap^/2- HcosA9 

A9 , for e > Ap^/2 - H cos A9 . 

(14) 

In the steady state, the {9,p) distribution is such that, 
for any given energy, all the microstates corresponding 
to that energy are equally probable. The boundaries on 
{9,p) imposed by the initial waterbag are no longer valid. 
Thus, the steady state distribution P{9,p) may be ex- 
pressed as 

where Qe{e{9,p)) is given by P^{e{9,p)), evaluated by 
Eq. (in]), with 6*1 = and 6*0 given in ^ with A9 re- 
placed by TT (see Fig.[T^). 
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FIG. 1: a) Graphical representation of the computation of 
Eq. (|12[) . The full lines represent the domain on which the 
integral is performed in order to compute for values of 
the energy e — 0, H and 2H. The boundaries 9q and 9i are 
displayed for each energy, b) Magnetization in the steady 
state nis vs. Ap for an initial waterbag state with A6 — 1. 
The full line is the theoretical prediction and the points are 
the results of numerical integrations of the Vlasov equation 
for the HMF model in Eqs. ((2} and ([3)|. Error bars quantify 
fluctuations of in the steady state. 



Integrating over p and using ([5]), it is straightforward 
to express, without any approximation, the marginal in 
6 as 



/oo 
de 
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The interpretation of is straightforward. For any 
given energy e, it is proportional to the ratio of the 
time taken to cover the orbit in the (0,p)-plane confined 
by the waterbag domain to the time corresponding to 
the orbit unrestricted by the initial condition. van- 
ishes for e < —H, becomes l/{AA9Ap) in the interval 
-H <€<€*, with e* = min(-i?cosA6', Ap2/2 - H), 
and monotonously decreases to zero in e* < e < — 
ifcosA6'. 



III. APPLICATION TO WATERBAG STATES 
OF THE HMF MODEL 

We now use these results in order to probe the steady 
states of the HMF model, where the particles are now 
interacting. For the distribution (|15p to be a stationary 
state of the corresponding HMF model one has to de- 
mand that the steady state magnetization of the model 
with non interacting particles satisfies H = m [isj , where 



dePg{e,H) cost 



(17) 



The strategy we propose is to solve the above self- 
consistency equation in H and then substitute the se- 
lected value of H into Eq. (jTS]) to get the steady state 
distributions, not only of the uncoupled system but also 



of the HMF. This is the main result of this paper. In 
applying these results to the HMF model, it is assumed 
that the magnetization takes its steady state value in- 
stantaneously a,t t — 0. This is of course not correct 
since the magnetization evolves towards its steady state 
value on a finite time scale, during which the momen- 
tum distribution also changes. Thus the correspondence 
between the initial state and the steady state distribu- 
tion is only approximate. This approximations works well 
for some initial conditions and not so well for others, as 
discussed below. 

We proceed by solving the self-consistency equation 
P?)) . For all A9 < tt and at small H the leading term in 
the marginal PeiO, H) behaves like y/H. This can be eas- 
ily proven numerically but also analytically by splitting 
the integral in (IT6l) in the three relevant domains of inte- 
gration [-H cos 9, -H cos AO], [-H cos A9, Ap'^ /2 - H], 
[Ap'^/2 - H, Ap'^/2 - iJ cos A9]. The integration on both 
the first and the second domain displays a \/H behavior 
and the first integral can be explicitly computed, giving 
2 yj2H{cos9 — cos A9). The third integral gives an order 
contribution. When computing m in formula ([TT]) the 
first and second integral contribute terms that are oppo- 
site in sign for large enough A9 and do not cancel each 
other. When A9 reaches tt the ^fH terms cancel, giving 
rise in this limit to an order H global contribution (we 
give below explicitly the result for A9 = tt). Because of 
the Vh behavior, the H = solution is always unstable. 
Therefore the magnetization is expected to be non zero 
at any value of Ap [ij, [l^. The theoretical magnetic 
curve obtained by and ([T7)) is compared with the 
numerical solution of the Vlasov equation for AO = 1 in 
Fig. [T]d. The agreement is quite good, considering that 
the theoretical prediction has no free parameters. 

In Fig. [2] we compare the marginals in 9 and in p of the 
steady state distribution (|15p with those obtained numer- 
ically by direct integration of the Vlasov equation. Apart 
from the oscillations observed in the tails, which are due 
to the filamentation of the initial waterbag distribution, 
the agreement is quite good. 

Let us now consider initial homogeneous waterbag 
states, for which A9 — tt. In analysing Eq. ((TT)) we find 
that m = H = is a solution for any Ap. At low Ap 
there exists one additional solution with m ^ 0, while 
in a rather narrow intermediate range Ap* < Ap < Ap 
two solutions with nonzero magnetization are present. In 
Fig. [3^ we plot m{H) — H as a function of H for Ap in 
the intermediate region where one gets three solutions of 
the self-consistency equation, two stable (S) and one un- 
stable (U). In Fig.[3)D we plot the non zero self-consistent 
steady value of the magnetization, ms, versus Ap. In 
this case the theoretical prediction does not agree with 
the numerical simulations of the Vlasov equation. We 
argue that this is due to the fact that, during time evolu- 
tion, there is a significant mixing of energy levels, due to 
strong variations of the effective field near a discontinu- 
ous first order transition [17|. 
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FIG. 2: Marginals in 6, Pe{9) = J dp P{9,p), (panel a) and 
in p, Ppijp) = f d6 P{6,p), (panel b) of the steady state dis- 
tribution (fT5)) for A9 ^ 1 and Ap = 1. The full lines are 
the theoretical predictions and the points are the results of 
numerical integrations of the Vlasov equation for the HMF 
model in Eqs. Q and ©. 




H Ap 

FIG. 3: a) m{H) ~ H vs. H for Ap = 0.72 within the range 
of existence of three solutions l/\/2 = Ap* < Ap < Ap ~ 
0.735. Two stable fixed point (S) with negative slopes and 
an unstable one (U) with a positive slope, b) The non zero 
steady state magnetization ms as a function of Ap. The 
dotted segment corresponds to unstable solutions and the full 
line to stable ones. The dot on the real axis marks the stability 
limit of the homogeneous state Ap*. The points are results 
of numerical integrations of the Vlasov equation for the HMF 
model in Eqs. (O and ®. 



IV. STABILITY OF THE HOMOGENEOUS 
WATERBAG STATE 

We now analyze the stability of the stationary solu- 
tions. For the zero magnetization state the stability limit 
is determined by the slope of m{H) at = 0, which is 
obtained by taking the small H limit in the integral (IT51) . 
We therefore expand Pg{6,H) in powers of H keeping 
only the linear term in H which is proportional to cos 6. 
This term is found to be H cos9/{2ttAp^). It is obtained 
by carrying out the integration in the energy domain 
—HcosO < e < — H, where Pe(e) is constant. 

The stability threshold for the m = solution is thus 
Ap* = l/\/2, as displayed in FigOjD. For Ap > Ap* 
the homogeneous state is stable and it becomes unsta- 



ble below Ap*. Using this value in the HMF model we 
get the energy (Ap*)^/6 -I- 1/2 = 7/12, which coincides 
with the threshold energy calculated using other meth- 
ods 0, Q . This analysis can be generalized to a generic 
initial momentum distribution. 

The stable inhomomogeneous solutions found for dif- 
ferent values of Ap are exact stable stationary solutions 
of the Vlasov equation, but they evolve in time for fi- 
nite N . Numerical studies of the time evolution of the 
magnetization are summarized in Fig. |4l where m{t) is 
displayed for several values of N . We obtain a data col- 
lapse when time is rescaled by 1/iV, implying that the 
relaxation time increases as 

ATT with 7 = 1 dllll. This 
is at variance with the relaxation of homogeneous states 
where 7 ~ 1.7 H US [2l|. 
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FIG. 4: Magnetization vs. time of the HMF model starting 
with the stable inhomogeneous steady state at Ap = 0.7 for 
iV = 2"(100), 2"(100), 2^2(50), 2^3(50), 2^''(20), 2^^{2Q) from 
left to right (in parenthesis the number of realisations). In 
the inset the same data are plotted as a function of t/N . 



V. CONCLUSIONS 

In summary, our approach allows one to derive exact 
inhomogeneous steady states. It gives approximate cor- 
respondence between initial and final states, which, for 
the HMF model, yields good results for initial inhomo- 
geneous waterbags. Moreover, one obtains the exact sta- 
bility limit of the zero magnetization state, since the wa- 
terbag initial condition is in itself a steady state with 
zero magnetization and the dynamics of each particle is 
therefore exactly decribed by the independent particles 
model dSl). 



VI. ACKNOWLEDGEMENTS 

We thank O. Cohen, O. Hirschberg, S. Gupta and V. 
Rittenberg for helpful discussions. The support of the 
Israel Science Foundation (ISF), the Minerva Founda- 
tion, with funding from the Federal German Ministry for 
Education and Reseach, INFN and the Belgian federal 



5 



government (lAP Project "NOSY" P6/02) is gratefully Seience and is also part of the ANR-lO-CEXC-010-01, 
acknowledged. This work was carried out while S.R. was Chaire d'Excellence project. 
Weston Visiting Professor at the Weizmann Institute of 



[1] T. Dauxois, S. Ruffo and L. Cugliandolo (Eds.), Long- 
Range Interacting Systems (Oxford University Press, [10 
(2009); F. Bouchet, S. Gupta and D. Mukamel, Phys- 
ica A 389, 4389 (2010); A. Campa, T. Dauxois and S. [11 
RuflFo, Phys. Rep. 480, 57 (2009). 

[2] Y. Y. Yamaguchi, J. Barrc, F. Bouchet, T. Dauxois and [12 
S. Ruffo, Physica A 337, 36 (2004). 

[3] M. Antoni and S. RufFo, Phys. Rev. E 52, 2361 (1995);P. [13 
H. Chavanis, J. Vatteville and F. Bouchet, European 
Physical Journal B46, 61 (2005). [14 

[4] D. R. Nicholson, Introduction to Plasma Theory, (John 
Wiley New York, 1983). 

[5] K. Jain, F. Bouchet and D. Mukamel, J. Stat. Mech. 11, [15 
8 (2010); J. Barre, A. Olivetti and Y. Y. Yamaguchi, J. [16 
Stat. Mech. P08002 (2010); A. Campa and P.-H. Cha- [17^ 
vanis, J. Stat. Mech. P06001 (2010); R. Bachclard, F. 
Staniscia, T. Dauxois, G. De Ninno and S. Ruffo, J. Stat. [18 
Mech., P03022 (2011). [19 

[6] I. Bernstein, J. M. Greene, M. D. Kruskal Phys. Rev., 

108, 546 (1957). [20 

[7] Y. Y. Yamaguchi, Phys. Rev. E, 84, 016211 (2011). 

[8] Y. Pomeau, Statistical mechanics of gravitational plas- [21 
mas, Lecture Notes, 2nd Warsaw School of Statistical 
Physics (2007). 

[9] X. Leoncini, T. Van Den Berg and D. Fanelli, Europhys. 



Lett. 86, 20002 (2009). 

P. de Buyl, D. Mukamel and S. Ruffo, Phil. Trans. R. 
Soc. A, 369, 439 (2011). 

D. Lynden-BeU, Mon. Not. R. Astron. Soc. 136, 101 

(1967). 

J. Barrc, T. Dauxois, G. De Ninno, D. Fanelli and S. 
Ruffo, Phys. Rev. E 69, 045501 (2004). 
A. Antoniazzi, D. Fanelli, S. Ruffo, and Y. Y. Yamaguchi, 
Phys. Rev. Lett. 99, 040601 (2007). 

Y. Levin, R. Pakter and T. N. Teles, Phys. Rev. Lett. 

100, 040604 (2008); T.N. Teles, Y. Levin, R. Pakter and 

F.B. Rizzato, J. Stat. Mech. P05007 (2010). 

F. Bouchet, Physica D, 237, 1976 (2008). 

M.-C. Firpo, Europhys. Lett. 88, 30010 (2009). 

R. Pakter and Y. Levin, Phys. Rev. Lett. 106, 200603 

(2011). 

P.-H. Chavanis J. Stat. Mech., P05019 (2010). 

M Joyce, T Worrakitpoonpon J. Stat. Mech, P10012 

(2010). 

L. G. Moyano and C. Anteneodo, Phys. Rev. E 74, 
021118 (2006). 

S. Gupta and D. Mukamel, Phys. Rev. Lett. 105, 040602 
(2010). 



